program nc_reducer

    use iasing_osse
    use netcdf

    implicit none
    ! Program counters 
    integer :: iowrite
    integer :: jo, jl
    integer :: nobs
    integer, allocatable :: listObs(:)

    ! Netcdf stuff
    character(len=64) :: ncfile 
    integer :: ncid, ouid
    
    ! Interpolation
    !real(JPR4) :: vprof(51)
    real(JPR8) :: dpaux

    ! iasing_osse stuff
    type(r4a_ids) :: db_ids, ou_ids
    type(r60_dat) :: db_dat
    type(r51_dat) :: ou_dat

    ! Interfaces
    include "checknc.interface"
    include "r60_initializer.interface"
    include "r60_reader.interface"
    include "r51_writer.interface"
    include "r51_idscreator.interface"
    include "r60tor51.interface"
    include "rad2brt_f.interface"

    ! Read selected
    open(1,file='selected',status='old',action='read')
    read(1,'(A)') ncfile
    read(1,*) nobs
    allocate(listObs(nobs))
    read(1,*) listObs
    close (1)

    ! Create new file
    call checknc(nf90_create("example_51.nc", NF90_HDF5, ouid))
    ou_ids = r51_idscreator(ouid, nobs)
    call checknc(nf90_enddef(ouid))

    ! Print info
    print*, ncfile

    ! Open old netcdf Databank
    call checknc(nf90_open(ncfile, NF90_NOWRITE, ncid))
    db_ids = r60_initializer(ncid)

    do jo = 1, nobs
        ! Read data
        db_dat  = r60_reader(ncid, db_ids, listObs(jo))
        
        ! Copy to r51_dat structure
        ou_dat = r60tor51(db_dat)

        ! Change units of gases vertical profiles and interpol
        do jl = 1, NLM51
            call rttov_q2v(gas_unit_specconc, real(ou_dat%vpr%q(jl),8), &
                 gas_id_co, real(ou_dat%vpr%co(jl),8), dpaux)
            ou_dat%vpr%co(jl) = dpaux
            call rttov_q2v(gas_unit_specconc, real(ou_dat%vpr%q(jl),8), &
                 gas_id_ch4, real(ou_dat%vpr%ch4(jl),8), dpaux)
            ou_dat%vpr%ch4(jl) = dpaux
            call rttov_q2v(gas_unit_specconc, real(ou_dat%vpr%q(jl),8), &
                 gas_id_co2, real(ou_dat%vpr%co2(jl),8), dpaux)
            ou_dat%vpr%co2(jl) = dpaux
            call rttov_q2v(gas_unit_specconc, real(ou_dat%vpr%q(jl),8), &
                 gas_id_ozone, real(ou_dat%vpr%o3(jl),8), dpaux)
            ou_dat%vpr%o3(jl) = dpaux
            call rttov_q2v(gas_unit_specconc, real(ou_dat%vpr%q(jl),8), &
                 gas_id_so2, real(ou_dat%vpr%so2(jl),8), dpaux)
            ou_dat%vpr%so2(jl) = dpaux
            call rttov_q2v(gas_unit_specconc, real(ou_dat%vpr%q(jl),8), &
                 gas_id_watervapour, real(ou_dat%vpr%q(jl),8), dpaux)
            ou_dat%vpr%q(jl) = dpaux
        end do

        !if (any(ou_dat%vpr%tem>TMAX) .or. any((ou_dat%vpr%tem<TMIN))) &
        !    print*, "Error in Temperature profile. Local index ", jo, listObs(jo)
        !if (any(ou_dat%vpr%q>QMAX) .or. any((ou_dat%vpr%q<QMIN))) &
        !    print*, "Error in Humidity profile. Local index ", jo, listObs(jo)
        !if (any(ou_dat%vpr%o3>O3MAX) .or. any((ou_dat%vpr%o3<O3MIN))) &
        !    print*, "Error in Ozone profile. Local index ", jo, listObs(jo)
        !if (any(ou_dat%vpr%co2>CO2MAX) .or. any((ou_dat%vpr%co2<CO2MIN))) &
        !    print*, "Error in CO2 profile. Local index ", jo, listObs(jo)
        !if (any(ou_dat%vpr%co>COMAX) .or. any((ou_dat%vpr%co<COMIN))) &
        !    print*, "Error in CO profile. Local index ", jo, listObs(jo)
        !if (any(ou_dat%vpr%ch4>CH4MAX) .or. any((ou_dat%vpr%ch4<CH4MIN))) &
        !    print*, "Error in CH4 profile. Local index ", jo, listObs(jo)
        
        ! Change radiance units
        ! IASI
        if (any(ou_dat%rbt%ria<=0. .or. ou_dat%rbt%ria>1000.)) then
            ou_dat%rbt%ria = -999.9
        else if (any(ou_dat%rbt%ria<100)) then
            ou_dat%rbt%ria = rad2brt_f(NIASI, ou_dat%dim%wia, &
                ou_dat%rbt%ria)
        endif
        ! IASI-NG A
        if (any(ou_dat%rbt%rng<=0. .or. ou_dat%rbt%rng>1000.)) then
            ou_dat%rbt%rng = -999.9
        else if (any(ou_dat%rbt%rng<100)) then
            ou_dat%rbt%rng = rad2brt_f(NIANG, ou_dat%dim%wng, &
                ou_dat%rbt%rng)
        endif
        ! IASI-NG B
        if (any(ou_dat%rbt%ngb<=0. .or. ou_dat%rbt%ngb>1000.)) then
            ou_dat%rbt%ngb = -999.9
        else if (any(ou_dat%rbt%ngb<100)) then
            ou_dat%rbt%ngb = rad2brt_f(NIANG, ou_dat%dim%wng, &
                ou_dat%rbt%ngb)
        endif

        ! Write data
        iowrite = r51_writer(ouid, ou_ids, jo, db_ids, ou_dat)
        
        if (nobs>10.and.mod(jo,nobs/100).eq.0) print*, "Observation ", jo
    enddo

    ! Close Netcdf Databank
    call checknc(nf90_close(ncid))
    call checknc(nf90_close(ouid))

end program
